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Abstract. This paper develops a fully discrete modified characteristic finite element method for a coupled system 
consisting of the fully nonlinear Monge- Ampere equation and a transport equation. The system is the Eulerian formulation 
in the dual space for the B. J. Hoskins' semigeostrophic flow equations, which are widely used in meteorology to model 
slowly varying flows constrained by rotation and stratification. To overcome the difficulty caused by the strong nonlinearity, 
we first formulate (at the differential level) a vanishing moment approximation of the semigeostrophic flow equations, a 
methodology recently proposed by the authors |17l 118) , which involves approximating the fully nonlinear Monge- Ampere 
equation by a family of fourth order quasilinear equations. We then construct a fully discrete modified characteristic finite 
element method for the regularized problem. It is shown that under certain mesh and time stepping constraints, the 
proposed numerical method converges with an optimal order rate of convergence. In particular, the obtained error bounds 
show explicit dependence on the regularization parameter e. Numerical tests are also presented to validate the theoretical 
results and to gauge the efficiency of the proposed fully discrete modified characteristic finite element method. 
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1. Introduction. The semigeostrophic flow equations, which were derived by B. J. 
Hoskins [22], is used in meteorology to model slowly varying flows constrained by rotation 
and stratiflcation. They can be considered as an approximation of the Euler equations 
and are thought to be an efficient model to describe front formation (cf. [231 [ID])- Under 
certain assumptions and in some appropriately chosen curve coordinates (called 'dual 
space', see Section |2|, they can be formulated as the following coupled system consisting 
of the fully nonlinear Monge-Ampere equation and the transport equation: 

(1.1) det{D^ip*) = a in x (0, T], 

dot 

(1.2) — + div(va) = inM^x(0,r], 

(1.3) a{x, 0) = ao in x {t = 0}, 

(1.4) Vil)* C fi, 

and 

(1.5) V = (v^* - x)^ = (v-:, - - 0). 

Here, f2 C is a bounded domain, a is the density of a probability measure on M^, and 
-0* denotes the Legendre transform of a convex function ip. For any w = {wi,W2-,w^), 
w"'- := (w2,— Wi,0). We note that none of the variables a,-?/^*, and v in the system is 
an original primitive variable appearing in the Euler equations. However, all primitive 
variables can be conveniently recovered from these non-physical variables (see Section |2] 
for the details). 



In this paper, our goal is to numerically approximate the solution of (1.1)-(1.5). By 
inspecting the above system, one easily observes that there are three clear difficulties for 
achieving the goal. First, the equations are posed over an unbounded domain, which makes 
numerically solving the system infeasible. Second, the ?/'*-equation is the fully nonlinear 
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Monge-Ampere equation. Numerically, little progress has been made in approximating 
second order fully nonlinear PDEs such as the Monge-Ampere equation. Third, equation 
(1.4) imposes a nonstandard constraint on the solution ip* , which often is called the second 
kind boundary condition for ip* in the PDE community (cf. [3l ITO]). 



As a first step to approximate the solution of the above system, we must solve (|l.l|) 
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(1.3) over a finite domain, U C M , which then calls for the use of artificial boundary 
condition techniques. For the second difficulty, we recall that a main obstacle is the fact 
that weak solutions (called viscosity solutions) for second order nonlinear PDEs are non- 
variational. This poses a daunting challenge for Galerkin type numerical methods such 
as finite element, spectral element, and discontinuous Galerkin methods, which are all 
based on variational formulations of PDEs. To overcome the above difficulty, recently we 
introduced a new approach in [T71 UHl [191 ISOl l25] , called the vanishing moment method 
in order to approximate viscosity solutions of fully nonlinear second order PDEs. This 
approach gives rise a new notion of weak solutions, called moment solutions, for fully 
nonlinear second order PDEs. Furthermore, the vanishing moment method is constructive, 
so practical and convergent numerical methods can be developed based on the approach 
for computing viscosity solutions of fully nonlinear second order PDEs. The main idea of 
the vanishing moment method is to approximate a fully nonlinear second order PDE by a 
quasilinear higher order PDE. In this paper, we apply the methodology of the vanishing 



moment method, and approximate (1.1)-(1.3) by the following fourth order quasi-linear 



system: 








(1.6) 




= a' 


in f/ X (0,r], 


(1.7) 






in f/ X (0,T], 


— + div (v"a") 


= 


(1.8) 


a^(a;,0) 


= ao{x) 


in X {t = 0}, 


where 








(1.9) 


:= {V^' - x)^ 


= (C 


- X2,Xi-^l^,0). 



It is easy to see that (1.6)-(1.9) is underdetermined, so extra constraints are required in 
order to ensure uniqueness. To this end, we impose the following boundary conditions 
and constraint to the above system: 



(1.10) 
(1.11) 
(1.12) 



du 

du 

ip^dx 

u 











on dU X (0,T], 
on dU X (0,T], 
t G (0,T], 



where v denotes the unit outward normal to dU . We remark that the choice of (1.11) 
intends to minimize the boundary layer due to the introduction of the singular perturba- 
tion term in (1.6) (see [17] for more discussions). Boundary condition (1.10) is used to 



minimize the "reflection" due to the introduction of the flnite computational domain U . 
It can be regarded as a simple radiation boundary condition. An additional consequence 
of (1.10) is that it also effec tively overcomes the third diffic ulty, which is ca used by the 
nonstandard constraint (1.4), for solving system (1.1)-(1.5). Clearly, (1.12) is purely a 



mathematical technique for selecting a unique function from a class of functions differing 
from each other by an additive constant. 



A characteristic FEM for the semigeostrophic flow 



3 



The specific goal of tliis paper is to f ormu late and analyze a modified characteristic 
finite element method for problem (1.6)-(1.12). The proposed method approximates the 
elliptic equation for ip^ by conforming finite element methods (cf. [8J) and discretizes 
the transport equation for by a modified characteristic method due to Douglas and 
Russell pLSj . We are particularly interested in obtaining error estimates that show explicit 
dependence on e for the proposed numerical method. 

The remainder of this paper is organized as follows. In Section [2| we introduce the 
semigeostrophic flow equations and show how they can be formulated as the Monge- 
Ampere/transport system (1.1)-(1.5). In Sectio n [3| we apply the methodology of the 
vanishing moment method to approximate (1. !)-(!. 5) via (1.6)-(1.12), prove some prop- 



erties of this approximation, and also state certain assumptions about this approximation. 
We then formulate our modified characteristic finite element method to numerically com- 
pute the solution of (|1.6)-(1.12). Section |4| mirrors the analysis found in [2Dj where we 



analyze the numerical solution of the Monge- Ampere equation under small perturbations 
of the data. Section |4] is of independent interests in itself, but the main results will prove 
to be crucial in the next section. In Section [5} under certain mesh and time stepping 
constraints, we establish optimal order error estimates for the proposed modified charac- 
teristic finite element method. The main idea of the proof is to use the results of Section 
|4]and an inductive argument. Finally, in Section [6} we provide numerical tests to validate 
the theoretical results of the paper. 

Standard space notation is adopted in this paper, we refer to [U EH [8] for their exact 
definitions. In particular, (-,■) and (■,■) denote the L^-inner products on U and dU , 
respectively. C is used to denote a generic positive constant which is independent of e 
and mesh parameters h and At. 

2. Derivation of the Monge- Ampere/transport formulation for the semigeostrophic flow 
equations. For the reader's convenience and to provide necessary background, we shall 
first give a concise derivation of the Hoskins' semigeostrophic flow equations ^22j and then 
explain how the Hoskins' model is reformulated as a coupled Monge- Ampere/transport 
system. Although our derivation essentially follows those of [22| ITOt [3] , we shall make an 
effort to streamline the ideas and key steps in a way which we thought should be more 
accessible to the numerical analysis community. 

Let f2 C denote a bounded domain of the troposphere in the atmosphere. It is 
well known [24j that if fluids are assumed to be incompressible, their dynamics in such a 
domain f2 are governed by the following incompressible Boussinesq equations which are a 
version of the incompressible Euler equations: 



(2.1) 

(2.2) 

(2.3) 
(2.4) 



Dn 



e 



DO 

In 

divu 






u = 



ge-i in(^x(0,T], 

in X (0,T], 

in X (0,T], 
on dn X (0,T], 



where 63 := (0, 0, 1), u = (ui, M2, M3) is the velocity field, p is the pressure, 6 either denotes 
the temperature (in the case of atmosphere) or the density (in the case of ocean) of the 
fluid in question. 9q is a reference value of 9. Also 
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denotes the material derivative. Recall that := ('U2,— mi,0). Finally, /, assumed to 
be a positive constant, is known as the Coriolis parameter, and g is the gravitational 
acceleration constant. We note that the term /u-*- is the so-called Coriolis force which is 
an artifact of the earth's rotation (cf. 130)'). 



Ignoring the (low order) material derivative term in (2.1) we get 



(2.5) 
(2.6) 

where 



dp 9 



9, 



V 



H '■- 



d d 



dxi ' 8x2 ' 



Equation (2.5) is known as the geostrophic balance, which describes the balance between 



the pressure gradient force and the Coriolis force in the horizontal directions. Equation 
(2.6) is known as the hydrostatic balance in the literature, which describes the balance 



between the pressure gradient force and the gravitational force in the vertical direction. 
Define 



(2.7) 



% := /"'(Vp)^ and 



ag 



U 



U 



3' 



which are often called the geostrophic wind and ageostrophic wind, respectively. 

The geostrophic and hydrostatic balances give very simple relations between the pres- 
sure field and the velocity field. However, the dynamics of the fluids are missing in the 
description. To overcome this limitation, J. B. Hoskins [22] proposed so-called semi- 



geostrophic approximation which is based on replacing the material derivative term 



by in ( |2.1 ). This then leads to the following semigeostrophic flow equations (in the 
primitive variables): 



Du 

Dt 



(2.8) 
(2.9) 

(2.10) 

(2.11) 
(2.12) 



Dt 



+ (Vp)^ = /u^ 
dp 9 



dxs 
DO 

Dt 
divu 

u 










in X (0,T], 
inVlx (0,T], 

in n X (0,T], 

innx (0,T], 
on dn X (0,T]. 



It is easy to see that after substituting Ug = f ^(Vp)"*", ( 2.8[ ) is an evolution equation 



for (Vp)"*". There are no explicit dynamic equations for u in the above semigeostrophic 
flow model. Also, by the definition of the material derivative, = ^ + (u • V)ug. 
We note that the full velocity u appears in the last term. Should u ■ V be replaced by 
Up ■ V in the material derivative, the resulting model is known as the quasi- geostrophic 
flow equations (cf. |24j ) . 

Due to the peculiar structure of the semigeostrophic flow equations, it is difficult to 
analyze and to numerically solve the equations. The first successful analytical approach 
is the one based on the fully nonlinear reformulation (1.1)-(1.5), which was first proposed 
in [5j and was further developed in [31 [23j (see ^HJ for a different approach). The main 
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idea of the reformulation is to use time-dependent curved coordinates so the resulting 
system becomes partially decoupled. Apparently, the trade-off is the presence of stronger 
nonlinearity in the new formulation. 

The derivation of the fully nonlinear reformulation (1.1)-(1.5) starts with introducing 
the so-called geopotential and geostrophic transformation 

(2.13) To "I — ki^p! $ := V'V'; where Xj^ := (xi, X2, 0). 



A direct calculation verifies that 



^ ■=xh + -pi^p)- 
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e 



Oof 



consequently, (2.8)-(2.10) can be rewritten compactly as 



(2.14) 
where 



Dt 



fJ{^-x) 



J 




For any a; G fi, let X{x,t) denote the fluid particle trajectory originating from x, i.e., 

^^^^ = u(X(x,t),t) Vt>0, 
X(x,0) = X. 
Define the composite function 

(2.15) vl/(.,t) := $(.,t) oX(-,t) = <f(X(-,t),t) = V^(X(-,t),t). 



Then we have from (|2.14|) 
(2.16) 



a*(x,t) 



dt 



fj{^{x, t) - X(x, t)) = fi^ix, t) - X{x, t))- 



Since the incompressibility assumption implies X is volume preserving, 

det(VX) = 1, 

which is equivalent to 



(2.17) 



g{X{x,t))dx= / g{x)dx G (7(fi). 



To summarize, we have reduced (2.8)-(2.11) into (2.15)-(2.17). It is easy to see that 
"9{x,t) is not unique because one has a freedom in choosing the geopotential ip- However, 
CuUen, Norbury, and Purser [12] (also see [TOl IH |23]) discovered the so-called Cullen- 
Norbury- Purser principle which says that \Ef(x, t) must minimize the geostrophic energy 
at each time t. A consequence of this minimum energy principle is that the geopotential 
ijj must be a convex function. Using the assumption that ip is convex and Brenier's polar 
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factorization theorem [5], Brenier and Benamou [3] proved e xisten ce of s uch a convex 
function ip and a measure preserving mapping X which solves (2.15)-(2.17). 



To relate (2.15)-(2.17) with (1.1), (1.2), and (1.4), let a{y,t)dy be the image measure 
of the Lebesgue measure dx by "^{x, t), that is 



g{^!{x,t))dx 



We note that the image measure a{y, t)dy is the push-forward \E'#(ix of dx by \E'(x, t), and 
a{y,t) is the density of "^^dx with respect to the Lebesg ue m easure dy. 

Assume that ip is sufficiently regular, it follows from (2.15) and (2.17) that 



(2.18) / g{^{x,t))dx= / g{Vij{X{x,t),t))dx= / g{Vij{x,t))dx V^? G Ce(M'). 

Using a change of variable y = Vipixjt) on the right and the definition of a{y,t)dy on 
the left we get 



g{y)a{y,t)dy 



giy)det{D'r{y,t))dy \/geC, 



where ip* denotes the Legendre transform of ip, that is, 
(2.19) 



Hence 



^* (y, t) = sup (x ■ ?/ - i){x, t)) . 



a{y,t) = dei{D'r{y.t)) 



which yields (1.1). 



For convex function ip^ by a property of the Legendre transform we have 'Vip*{y,t) = 
X eQ. Hence Vip* C Q, therefore, ( L4|) holds. 

Finally, for any w G C^([— 1, T]; M'^), it follows from integration by parts and (2.16) 
that 



w(^(x,0),0) dx 
n Jo Jn 

T 



dw{^!{x,t),t) 
Jt 



dxdt 



^ , , d-^ix^t) dw(^(x,t),t) , , , 
Vw{^{x, t),t) + ^ i ' ^' ^ !> dxdt 



Jn 

T 



dt 



dt 



[ [ \Vw{^{x,t),t)-f{^{x,t)-X{x,t))^ + 
Jo Jn ^ 



dw{'^{x,t),t) 



dxdt. 



Making a change of variable y = Vip{x, t) and using the definition of a{y, t)dy we get 



(2.20) 



Jm? 



f dw{y,t) 
,1 dt 



+ f^{y,t) ■Vw{y,t)\a{y,t) dydt + / w{y,0)a{y,0) dy = 0, 
^ Jm.3 



where v is as in (1.5). Hence, 

dt 



+ /div(v(|/,t)«(y,t)) =0, 



which gives (1.3) as / = 1 is assumed in Section [T] 
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We remark that (2.18) and (2.20) are weak formulations of (1.1) and (1.2), respectively. 



We also cite the following existence and regularity results for (1.1)-(1.3) and refer the 
reader to [3] for their proofs. 

Theorem 2.1. Let Qo,Q C be two bounded Lipschitz domain. Suppose further 
that ao G Lp(M^) with oq > 0, supp{ao) C VLq, and J^^ao{x)dx = Then for any 



T > 0, p > 1, ( 1.1[ )-(1.3) has a weak solution (■?/'*,«) in the sense of (2.18) and ( |2.20 ). 
Furthermore, there exists an R> such that supp{a{x,t)) C -B_r(0) for all t G [0,TJ and 



aeL^{[0,T];LP{Bnm) 

G L°°([0,T];iyi'°°(fi)) 
i/j* G L°°([0,T];iyi'°°(M3) 



nonnegative, 

convex in physical space, 

convex in dual space. 



Remark 2.1. (a). The above compact support result for a justifies our approach of 
solving the original infinite domain problem on a truncated computational domain U , in 
particular, if U is chosen large enough so that -B_r(0) C U. 

(b). Since a and ip* are not physical variables, one needs to recover the physical 
variables u and p from a and ip* ■ This can be done by the following procedure. First, 
one constructs the geopotential ip from its Legendre transform ip* . Numerically, this can 
be done by fast inverse Legendre transform algorithms. Second, one recovers the pressure 
field p from the geopotential ip using (2.13). Third, one obtains the geostrophic wind \ig 



and the full velocity field u from the pressure field p using (2.7). 

(c) . Recently, Loeper J2^ generalized the above results to the case where a is a global 
weak probability measure solution of the semigeostrophic equations. 

(d) . As a comparison, we recall that two-dimensional incompressible Euler equations 
(in the vorticity-stream function formulation) has the form 



A(p = uj znfix(0,T], 
+ div(ucu)=0 znfix(0,T], 



u 



Clearly, the main difference is that (p-equation above is a linear equation while ip* in (1.2) 
is a fully nonlinear equation. 

We conclude this section by remarking that in the case that the gravity is omitted, 
then the flow becomes two-dimensional. Repeating the derivation of this section and 
dropping the third component of all vectors, we then obtained a 2-d semigeostrophic flow 
model which has exactly the same form as (1.1)-(1.5) except that the definition of the 
operator (■)-'- becomes w-*- := {w2, —Wi) for w = (^1,^2), and v in (1.5) is replaced by 



(^2 -X2,Xi - 



Similarly, v*^ in (1.9) should be replaced by 



In the remaining of this paper we shall consider numerical approximations of both 2-d 
and 3-d models. 



3. Formulation of the numerical method. 
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3.1. Formulation of the vanishing moment approximation. As pointed out in Section [T| 
the primary difficulty for analyzing and numerically approximating the semigeostrophic 
equations (1.1)-(1.5) is caused by the strong nonlinearity and non-uniqueness of the ip*- 
equation (i.e., Monge-Ampere equation, cf. [HEl])- The strong nonlinearity makes the 
equation non-variational, so any Galerkin type numerical methods is not directly appli- 
cable to the fully nonlinear equation. Non-uniqueness is difficult to deal at the discrete 
level because no effective selection criterion is known in the literature which guarantees 
picking up the physical solution (i.e., the convex solution). Because of the above difficul- 
ties, very little progress was made in the past on developing numerical methods for the 
Monge-Ampere equation and other fully nonlinear second order PDEs (cf. [131 EHl 129]). 

Very recently, we have developed a new approach, called the vanishing moment method, 
for solving the Monge-Ampere equation and other fully nonlinear second order PDEs (cf. 
[TT] [T8| [19] [20] [25] [26]). Our basic idea is to approximate a fully nonlinear second order 
PDE by a singularly perturbed quasilinear fourth order PDE. In the case of the Monge- 
Ampere equation, we approximate the fully nonlinear second order equation 



(3.1) 



det{D^w) 



by the following fourth order quasilinear PDE 

-eA'^w' + det{D^w') 



{6>0) 



accompanied by appropriate boundary conditions. Numerics of fT8| [T9] [20] [25] show that 
for fixed > 0, converges to the unique convex solution w of (3.1 ) as e — *• 0"*". Rigorous 
proof of the convergence in some special cases was carried out in [17J. Upon establishing 
the convergence of the vanishing moment method, one can use various well-established 
numerical methods (such as finite element, finite difference, spectral and discontinuous 
Galerkin methods) to solve the perturbed quasilinear fourth order PDE. Remarkably, our 
experiences so far suggest that the vanishing moment method always converges to the 
physical solution. The success motivates us to apply the vanishing moment methodology 



to the semigeostrophic model (1.1)-(1.5), which leads us to studying problem (1.6)-(1.12). 

Remark 3.1. Since a perturbation term is introduced in (1.6), it is also natural to 
introduce a "viscosity" term —eAa on the left-hand side of (1.7). We believe this should 



be another viable strategy and will further explore the idea and compare the anticipated 
new result with that of this paper. 



Since (1.6)-(1.7) is a quasilinear system, we can define weak solutions for problem 
(1.6)-( l7l2p in the usual way using integration by parts. 

Definition 3.1. A pair of functions [il)' , a^) e L'^{{0,Ty, H^{U))xL\{0,Ty, H^{U))n 



H^{{0,T)] L'^{U)) is called a weak solution to (1.6)-(1.12) if they satisfy the following in- 
tegral identities for almost every t G (0,T).- 

-e{Atlj',Av) + {det{D^tlj'),v) = 



(3.2) 

(3.3) 

(3.4) 
(3.5) 

here 



dt 



w 



+ (v" • Va^w) 



Vx G L\U), 



-- 

= («o,x) 
= 0, 

3 and = {ip^^ — X2, Xi — ip^J when d = 2, 



{ipl,^ — X2,Xi — ipl^ , 0) when d = 
and we have used the fact that divv^ = 0. 

For the continuation of the paper, we assume that there exists a unique solution to 



(1.6 )-( 1.12) such that il)'^{x,t) is convex, a'^{x,t) > 0, and supp a^(x, t) C -B_r(0) C U for 
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all t e [0,T]. We also assume G L\{0,T); H'{U)) (s > 3), a' G L\{0,T); Hf{U)) 
(p > 2), and that the following bounds hold (cf. [TTj) for almost all t G [0, T] 



(3.6) 
(3.7) 



1-j ■ 



Oie — ) (j = 1,2,3), 
= 0(£i-^) (j = l,2). 



e 



\a 



L°° - 
14/1,00 



where $^ = cof{D'^ip^) denotes the cofactor matrix of D'^ip^. 

As expected, the proof of the above assumptions is extensive and not easy. We do 
not intend to give a full proof in this paper. However, in the following we shall present 
a proof for a key assertion, that is, a^(x,t) > in f/ x [0,T] provided that ao{x) > 
in M'^((i = 2,3). Clearly, this assertion is important to ensure that ip^{-,t) is a convex 
function for all t G [0,T]. 

Proposition 3.2. Suppose {a^yip'^) is a regular solution of (1.6)-(1.12). Assume 
ao{x) >OinR'^{d = 2, 3), then a^{x, t) > in U X [0, T] . 

Proof. For any fixed (x,t) G f/ x (0,T], let {x, t ; s) denote the characteristic curve 
passing through {x,t) for the transport equation ( |1.7 ), that is. 



dX^ix, t: s) 



y^%X'ix,t;s),s) 



ds 

X{x, t; t) = X. 

Then the solution at (x, t) can be written as 

a'{x,t) = ao{X'ix,t;0)). 

Hence, a^{x,t) > for all (x, t) G f/ x [0,T]. The proof is complete. □ 

3.2. Formulation of modified characteristic finite element method. Let 7/j be a quasiuni- 
form triangulation or rectangular partition of U with mesh size h G (0, 1) and V'^ C H'^{U) 
denote a conforming finite element space (such as Argyris, Bell, Bogner-Fox-Schmit, and 
Hsieh-Clough-Tocher finite element spaces [8] when d = 2) consisting of piecewise poly- 
nomial functions of degree r (> 4) such that for any v G H'^{U) {s > 3) 



(3.8) 



inf \\v - VhWw < ^WvWh" 



j = 0, 1,2; £ = min{r + l,s}. 



Also let be a finite dimensional subspace of H'{U)) consisting of piecewise polynomials 
of degree k {> 1) associated with the mesh 7/j. 
Set 



(3.9) := G V\ ^ 

(3.10) Wl^ ■.= {whe 







dU 



dU 



0}, 



fl,vn 



It is easy to check that 

d_ 



d_ 
di 



V 



r : = 



1 



1 + V 



e\2 



'I + |v£|2 \dt 



d 



+ v^- V 



Hence, from (1.7) we have 
(3.11) 



dr 



1 + |vs|2 V dt 



0. 
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Here we have used the fact that divv^ = 0. 

For a fixed positive integer M, let At := ^ and '■= rn^t for m = 0, 1, 2, ■ ■ ■ , M. 
For any a; G t/, let x := a; — v^(x,t)At. It follows from the Taylor's formula that (cf. 

mm) 



(3.12) 



tm) tm) - t^^i) 



dr 



At 



0{At) 



for m = 1, 2, 



.M. 



Borrowing the ideas of [T^ [15] , we propose the following modified characteristic finite 
element method for problem (1.6)-(1.12): 



Algorithm 1: 

Step 1: Let a° be the finite element interpolation or the elliptic projection of a^. 



Step 2: For m = 0, 1, 2, . . . M, find {ip'^, G x such that 



(3.13) -e{Ai,^,Avh) + (det(D2^;r), ^/.) = ("JT,^/.) + {e\vh) Vt;, G Vq", 

(3.14) (C1)=0, 



(3.15) 
where 



■= (3;;,) 



X - Y]^At, 



In the case that is the continuous linear finite element space (i.e., k - 
the following lemma. 

Lemma 3.3. Let k = 1 in the definition ofW^, suppose that a^>0 in 
then the solution of Algorithm 1 satisfies > in U for all m >1. 

Proof. In the case = 1, (3.15) immediately implies that 



1), we have 



a 



m+l 



where {Pj} denote the nodal points of the mesh and Pj := Pj — v™At. Suppose that 
d^{Pj) > for all j. Since the basis functions of the linear element are nonnegative, then 



we have a™(Pj) > for all j. Hence, 



a 



m+l 



(Pj) > for all j. Therefore, the assertion 



follows from the induction argument. □ 

Remark 3.2. The positivity of for all m > 1 gives hope to verify the convexity of 
il)"^, which remains as an open problem (cf. [7^ [7P|. WOp. For high order finite elements 
(i.e., k >2), might take negative values for some m > although we shall show later 
that the deviation from zero must be very small. 



Let (■?/''^,q;^) be the solution of ( 1.6 )-( 1.12) and {ip'^,a'^) be the solution of (3.13) 



(3.15). In the subsequent sections we prove existence and uniqueness for (t/;™ 



^) and 

under certain mesh 

and time stepping constraints. To this end, we first study (|3.13|) independently, which 



provide optimal order error estimates for i})'^{tm)—i''h ^^"^ oi'^{tn 



motivates us to analyze finite element approximations of the Monge-Ampere equation 
with small perturbations of the data. Such an analysis enables us to bound the error 



V'^(^m) ~ i^h" terms of of the error Q;^(t^ 



We use similar techniques to those 



developed in [20] to carry out the analysis. With this result in hand, we use an inductive 
argument in Section [sjto get the desired error estimates for both ^^{tm)—^'^' «^(^m) — 



Oil 
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4. Finite element approximations of the Monge- Ampere equation with small pertur- 
bations. As mentioned above, analyzing the error ip^{tm) — V'/T motivates us to consider 
finite element approximations of the following auxiliary problem: for e > 0, 



(4.1) 


-6A\^ + det(DV) 


= ^ (> 0) 


in f/, 


(4.2) 


du 


= 


on dU, 


(4.3) 


dAu'P 
du 


= e 


on dU, 


(4.4) 


K, 1) 


= 0, 





whose weak formulation is defined as seeking -u"^ G H'^{Q) such that 

(4.5) -e{Au'^, Av) + {det{D\'^),v) = v) + {e\ v) 'iv G H'^iU), 

(4.6) [u^, 1) = 0. 

We note that the finite element approximation of a similar Monge- Ampere problem was 
constructed and analyzed in [20], where the Dirichlet boundary condition was considered 
and the right-hand side function ip is the same in the finite element scheme as in the PDE 



problem. In this section, we shall study the finite element approximation of (4.1)-(4.4) in 



which if is replaced by := (p + 5(p, where 5(p is some small perturbation of (p. Specifically, 



we analyze the following finite element approximation of (4.1)-(4.4): find G such 
that 

(4.7) -e{AulAvh) + {dei{D\l),v,,) = {^,Vh) + {e\vh) ivj, G 



As expected, we shall adapt the same ideas and techniques as those of [20] to analyze 
the above scheme. However, we shall omit some details if they are same as those of [20] 
but highlight the differences if they are significant, in particular, we shall trace how the 
error constants depend on e and 6ip. Also, since the analysis in 2-d and 3-d are essentially 
the same, we shall only present the detailed analysis of the three dimensional case and 
make comments about the two dimensional case when there is a meaningful difference. 



To analyze scheme (4.7), we first recall that (cf. [20]) the associated bilinear form of 
the linearization of the operator M^{u'^) := —eA'^u'^ + det{D'^u'^) at the solution u'^ is 
given by 

(4.8) B[v, w] := e{Av, Aw) + (<l>'^Vt;, Vw), 

where = coi{D'^u'^) denotes the cofactor matrix of D'^u'^ . 

Next, we define a linear operator T*^ : —>■ Vi such that for Wh G V^, T'^{wh) G 
is the solution of following problem: 

(4.9) B[wh - T^{wh),Vh] = e{Awh, Av^) - idet{D^Wh),Vh) 
It follows from [20^, Theorem 3.5] that T'^ is well-defined. Also, it is easy to see that 



any fixed point of is a solution to (4.7). We now show that if ||5(y9||i2 is sufficiently 
small, then indeed, T'^ has a unique fixed point in a neighborhood of m*^. To this end, we 
set 



,(p) := {vh G V;^ \\vh - Ihu'^Wh^ < p], 
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where Ihu'^ denotes the finite element interpolant of m'^ onto V^. 

Before we continue, we state a lemma concerning the divergence row property of 
cofactor matrices. A short proof can be found in [16]. 

Lemma 4.1. Given a vector-valued function w = (wi, W2, • • • , w„) : U — >■ M" . Assume 
w G [C^(t/)]". Then the cofactor matrix cof{Dw) of the gradient matrix Dw ofw satisfies 
the following row divergence-free property: 

n 

(4.10) div {coi{Dw))i = ^ d^^{coi{Dw))ij = for i = 1, 2, ■ ■ ■ , n, 

where (cof(Dw))j and {cof{Dw))ij denote respectively the ith row and the {i, j)-entry of 
cof(L'w). 

Throughout the rest of this section, we assume u"^ G if*, set £ = min{r + 1, s}, and 



assume the following bounds (compare to those of [20] and (3.6)): for j = 1,2,3 



(4.11) IIm'^II^, =0(£^), \\u'^\\w2,^ = Oie-'), ||$'^||ioo = 0(£-^) 
We then have the following results. 

Lemma 4.2. There exists a constant Ci{e) = 0{e^^) such that 

(4.12) - T^hunU^ < Ci{e){e~^h'-'\\u''yi + \M\h-^)- 



Proof. To ease notation set Sh = IhU"^ ~ ^'''(iftW^) and rj = IhU"^ — -u*^. Then for any 
Vh G Vq, we use the Mean Value Theorem to get 

+ {(p,Vh) + {e^,Vh) 

= e{Ari, Avh) + (T^ : D^u'' - hu^), vh) + (5<^, vh), 

where = cof (rD^^j^^^) + (1 - r)D'^u^) for r G [0, 1]. 
On noting that 

|Ti,| = |cof(rD2(4^'^) + (i_^)z^V),,| = \det{TD\hu'')l^ + {l-r)D\\)\, 

where D'^u'^ | . . denotes the resulting 2x2 matrix after deleting the i''^ row and j*^ column 
of D'^u'^, we obtain 

\mi,\<2 max {\TiD\hu^)ke + il-r)iD\'^)ki\Y 

<C max |(Z}V)hP < C||Z^V||ioo. 



Hence, from (4.11) it follows that ||T^||ioo = 0{6 ^). Thus, 

B[sh,Vh\ < £:||A?7||z,2||At;,,||i2 + Ce''^\\D'^r]\\L2\\vh\\L2 + \\5ip\\H-^\\vh\\m 

< C{e-^7]\\H2 + \\6^\\H-^)\\Vh\\H^. 

Finally, using the coercivity of B[-, ■] we get 

\\sh\\H^ < Ci{e){e-'h'-^\\u''\\H^ + \M\h-^). 
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The proof is complete. □ 

Lemma 4.3. There exists Hq > such that for h < Hq, there exists an p = p{h,e) 
such that for any Vh,Wh € B/j(p) there holds 

(4.13) ||T'^K)-T'^K)||h2 < ^-\\vH-Wh\\m. 



Proof. From the definitions of T'^(%) and T'f{wh) we get for any Zh G Vq 

= {^""{Vvh - Vwh), Vzh) + (det(Z}2^^) _ dei{D^Wh), z^) . 

Adding and subtracting det(i5^f^) and det(i5^w^), where and denote the stan- 
dard mollifications of Vh and Wh, respectively, yields 

B[T^{vn)-T^{wh),zu] 

= {^'^{trr,){Vvh - Vwh), Vzh) + (det(D2<) - det(D2^;:), z^) 

+ {det{D\h) - det(D2^;:), z^) + {dei{D^w''^) - dei{D^Wh),Zh) 

+ {det{D\h) - det{D''v''^),Zh) + {dei{D^w''^) - dei{D^Wh),Zh), 

where = cof{DX + t{D^w';^ - D^)) for r G [0, 1]. 
Using Lemma |4.1| and Sobolev's inequality we have 

(4.14) B[T^v,)-T^Wf,),z,] 

= (($^ - ^h){VvH - Vwh),VzH) + {^hiVvH - Vv';^),Vzh) 
+ (^h(V< - Vwh), Zh) + {^^i{D\h) - det(DX). ^h) 

+ \\wh - w^Jh^] + ||det(DV) - det(DX)IU2 
+ lldet(DV) -det(D2^0|U2}||z^||j^2. 

It follows from the Mean Value Theorem that 

- ^hUIl^ = II det{D\%) - det{DXl^ + r(DXL, - D^I^Ml^ 
= \\A^^ : {D'u% - {DXl, + r{D'wf^l^ - I?XL,)))IU^, 

where A^^' = cof(D2^^| + X{D'^v^\.. + t{D'^w^^\.. - DX\..))) G R^^^ ^ ^ [g, 1]. 
We bound HA^-' Hioo as follows: 

IIA^^IU^ = ||cof(DV|^^. + X{DXl^ + r(DXL, - DXI.MIl- 
= \\D\% + A(DXL, + r{D'w^,l^ - DXl^))h^ 
< C{e-' + h-lp + WD'v';^ - DV)||loo + ||/^'< - D'whWl^), 
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where we used the triangle inequahty followed by the inverse inequality and (4.11). Com- 
bining the above two inequalities we get 



Hence, 



(4.15) 11$^ - < C{e-' + h-lp + \\DX - DWWl- + - D'w.h^) 



Applying (4.15) to (4.14) and setting yU — > yield 
Using the coercivity of B[-, ■] we get 

(4.16) ||T^(i;,)-T'^(w,)lli/^ <C£-^(£-^ + /i-tp)(/i^-l^i'^||^.+p)||i;,-ti;,||^2. 



Finally, setting ho = O 
from (4.16) that 



h < ho, and p = 0(min{£:^, }), it then follows 



\\T^{vH)-T'^{wh)\\H^ < -\\vj,-w4h^ 



The proof is complete □ 

With the help of the above two lemmas, we are ready to state and prove our main 
results of this section. 

Theorem 4.1. Suppose WdifWH-^ = 0(min{e^, 5^/12 }). Then the re e xists an hi > 
such that for h < hi, there exists a unique solution uf^ G Vi solving (4.7). Furthermore, 
there holds the following error estimate: 



(4.17) 



- uIWh^ < C2{s){E-^h'-^\\u^H^ + \m\H-^) 



With C2{e) = 0{e-^). 

Proof. To show the first claim, we set 



hi = O \ min 



Hi 



-2|U,'P| 



Fix h < hi and set pi = 2Ci{e)[e "^h^ ^\\u^\\hi 
Cmin{e^, e/ia }. 



H-'i)- Then we have pi < 



Next, let Vh € B/j(pi). Using the triangle inequality and Lemmas 4.2 and 4.3 we get 



<Ci{e){e-^h'-^\\u'^\\H^ 



<P1 + Pl 

-2 2 



Pi- 
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Hence, T'-^{vh) G IBh(pi)- In addition, by (4.13) we know that T"^ is a contracting mapping 
in B/i(pi). Thus, the Brouwer's Fixed Point Theorem |21| guarantees that there exists a 
unique fixed point G Mh{pi) which is a solution to (4.7). 
Finally, using the triangle inequality we get 

□ 

Theorem 4.2. In addition to the hypothesis of Theorem 4jJ_. assume that the lin- 
earization of at u"^ (see (4.8 )j is H^-regular with the regularity constant Cs{s). Fur- 
thermore, assume that ||5</)||//-2 = 0{C'g^{e)e'^). Then there exists an h2 > such that 
for h < h2, there holds 



(4.18) 



< Cs{6)\C3{6)h'-'\\u^\\H, + {C^{e)h + l)\\S<f\\ 



where 6*3(5) = C2{e)e 2 and C4{e) = C2{s)e 2 
Proof Let 



and M^'^ denote a standard mollification of u^. We note that 



e'^ satisfies the following error equation: 



(4.19) £(Ae'^, Azh) + (det(D"u^) - det(D^M'^), Zh) + {5ip, zh) = 
Using (4.19), the Mean Value Theorem, and Lemma [4T we have 

(4.20) = £(Ae^ Az,) - (<I>V«'^ - n^), Vz,) + (5<^, z,) 

+ (det(DX)-det(DX'^),^,), 

where i> = cof(D2u^''' + t{D'^u'^ - D^ul'^)) for r G [0, 1]. 

Next, let G Vo n be the unique solution to the following problem: 

= (Ve'^,Vz) yzeVo. 

The regularity assumption implies that 

(4.21) \mH^<Csis)\\\/e^L^. 
We then have 

= £(Ae'^,A0) + (<l>'^V0,VeJ^) 

= £(Ae'^, A(0 - 40)) + ($^Ve^, V(0 - 40)) + £(Ae^ A(40)) 
+ ($'^Ve^ V(40)) - £(Ae^ A(40)) - (<i>V(M'^ - <'^), V(40)) 



|Ve'^||i2 



- (det(D2^^) _ det(D^<'^), 40) - {6^, m 



(4.22) 



= e{Ael A(0 - 40)) + ($'^Ve^ V(0 - 40)) 
+ (($^ - <l>)Ve^ V(40)) + miut^" - <), V(40)) 
+ {det{D\''^'^) - det(D2M^), 40) - {Sif, Ih<l>) 

< £||Ae^|U2||A(0 - 40)||l2 + C||$'^|U2||e^||H2||0 - 40II//2 



+ C||$^-l>|U2||Ve'^|U2||V(40)IU-+C^||$||L2||Mr-<||H2||40||//2 



+ lldet(DX) -det(DX'^)llL2||40||L2 + \\6^h\\H-4h<P\\H^ 
-h\\e^\\H2 + 11$^ - $||L2||Ve^|U2 + \\^L4ut'' - 



+ lldet(DM) -det(DX'^)|U. + 11011^3 
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We bound as H^"^ - $||l2 as follows: 

||(<l>^ - l>),,|U2 = llcof(DV),, - cof{{D'ul'^ + t{D\^ - D'ut'n)i,)\\L- 



II det(D\'^ 
||A'^' : {D\ 



2 ip\ 



dV\ 



iD\t% + riD''u 



<2||A'iz.-(pV-DX 



v'll I II n'2„,'P 



L2 



where A*^ = cof{D^u'^\.. + XiD^u'^'^l.. + TiD^u"^].. - ^ e [0, 1]. Notice that 



we have abused the notation A*-' by defining it differently in two proofs. 
To estimate ||A*^||ioo, we note that A*^' G R^''^. Thus for h < hi 



\\a''\\l^ < c(||dV|u- + ||dV - z^MlU- + II^M - d 



U 



h 



) 



< C(6-' + C2ie){e-'h'-l\\u^j,. + h-l \\6^h\\h-^) + ll/^M'^ - ^XlU-) 



<C{e-' + \\D'ur-D'unL^). 



where we have used the triangle inequality, the inverse inequality, and (4.11). Therefore, 
(4.23) 11$^ - $^1^2 < C(^e~^ + p'^r - D^l 



'■h\\L° 



Using (4.23) and setting /i ^ in (4.22) yield 



llVe'^lli^ <c\^e-h\\e^\\H2 + 

+ e-'C2ie){e''h'-''\\u^\\Hi + \\5^h\\H-^)\\^e^L'^ 



HA. 



It follows from (4.21) that 



\\Ve^L^<C,{e)[e-^^h\\el\\H^ + \\5v\\L2 



Set ho 



O 



Cs{e)\\u^\\jjt 

h < min{/ii, /12} 



1 

1-2 



. On noting that ||(5v9||//-2 < £{Cs{e)C2{s)) ^ we have for 



llVe'^IU^ <C,(£)(C2(£)£-t/i^-iM'^||^^. + (£-^C2(£)/i + 1)||5v^,||h-2). 



Thus, (4.18) follows from Poincare's inequality. The proof is complete. □ 

Remark 4.1. Let {ip'^^a^) he generated by A lgo rithm 1. If ||a;^(tm) — "/TIU^ = 
0(min{e^, e^/ia , (7~^(5)£^}) , then by Theorems 4-1 oind 4-2. for h < min{/ii, /i2}, there 
exists a unique ip"^ G solving (3.13), where 
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Furthermore, we have the following error bounds: 

(4.24) llV^^(t^) -d|^,2 < C2{e){e-'h'-'\\r{trn)\\m + \\a'{tm) - a^i^), 

(4.25) ||^^(t™) - < Csis){C,is)h'-'U'{tJ\\H^ 

+ {C,{e)h + l)\\a%t^)-anL^). 

Remark 4.2. In the two dimensional case, 



1 

1-2 

hi = 0{ „ , „ ) , h2 = 0' ^ 



3 



Il2([0,T];H«) / ' \Cs{£)\\i''^\\L'2{[0,T];H<^) 



1 

e-2 



and (4.25) holds with Cs{e) = C2{e)e^^ and C4{e) = C2{s)e^^. Furthermore, we only 



require ||a^(tm) — "/TiU^ = 0(mm{e'^,C^ ^{e)e}). 

5. Error analysis for Algorithm 1. In this section we shall derive error estimates for 
the solution of Algorithm 1. This will be done by using an inductive argument based on 
the error estimates of the previous section. Before stating our first main result of this 
section, we cite the well-known error estimate results for the elliptic projection of a{tm), 
which we denote by xT ^ • Let := a(-, t^) ~Xh'{')^ then there hold the following 

estimates for a{tm) — X/T?"^ ^ 1 (cf- [4, 8J): 

(5-1) ||'^||L2([0,r];L2) + /i||t^||L2([0,r];/fi) < C k^aW L'2{[0,T];H3), 

||^t|U2([o,r];L2) + h\\uJt\\L2([o^T];m) < C*^'' || || L2 ([0,T] ) , 



where j := min{A; + l,p}. As in Section |4| we set i = min{r + 1, s}. 

Theorem 5.1. There exists > such that for h < m.m{hi, h2, h^} there exists 
Ati > such that for At < min{Ati, h"^} 

(5.2) max ||a^(t„) -a™||i2 < C^ie)! At\\al^\\L2(io,T]xRS) 

0<m<M (. 

L2([0,T];H<) (, 

(5.3) max ||^^(t„) - ^P^Wh'^ < Cjie)! At\\aUL2(^[QT]xR^) 

0<m<M (. 

+ h^[\\'^''\\L2{[0,T];Hf) + ||L2([o,T];//i)] 

+ Ce{e)h'~^r\\LH[onm)}^ 

(5.4) max \\i;%tm) - ipTWm < C8is)\ At\\aULmo,T]xR^) 

0<m<M K 

+ [\W\\L-i{[0,T];Hi) + || «f ||L2([o,T];//i)] 

+ C^{e)h'~^\\r\\L-(\0,T\;H')]: 
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where C^{e) = 0{e-^), C,{e) = Cs{e)C^{e), Cj{e) = C2{e)C,{E), Cs{e) = Cs{e)C,{E), 



and Cs{b) is defined in Theorem 4-2 



Proof. We break the proof into five steps. 

Step 1: The proof is based on two induction hypotheses, where we assume for m 
0,1,--- ,k, 



(5.5) 
(5.6) 



\\D'^nL<^ = o{E-'). 



We first show that the claims of the theorem hold when k = 0. Let 



= O \ min 



(— 



HI 



ll«o|| 



HI 



' \Csie)\\ao\\ 



Hi 



From (5.1 ) we have for h < 

Wao-ahh^ <Chmao\\H^ < Cmm{e\e''h'^,C;\e)e''}. 



By Remark 4.1 



there exists solving (3.13). On noting that hi < C 



we have for h < minj/ii, h2, h^} 



< C(e-' + h~IC2{e){e~^h'-^rm\H^ + h'lKyA ) < Ce-\ 



The remaining four steps are devoted to show that the estimates hold for m = k + 1. 
Step 2: Let ^™ := a™ — x™- By ( |3.15 ) and (1.7), and a direct calculation we get 



(5.7) 



(r+' -r,r^') = (At<(wi) - («^(wi) -«ut™)),r+') 



+ ( 



m+1 rr"* £-"1+1 



h 5 



where := ^™(xh), a^(tm)) := a''{xh,tm), and cu™ := uj"^{xh)- 

We now estimate the right-hand side of ( |5.7 ). To bound the first term, we write 



Atal{x,tm+i) - [a''{x,tra+i) - a'ixh^trr,)] 

= At al{x, tm+i) - [a^(a;, t^+i) - a%x, tm)] + [a%Xh, tm) - a%x, tm)] ■ 

Using the identity 

Atal{x,trn+i) - [a'{x,tm+i) -a^(^,tm)] 

= / V|x(r)-a;|2 + (t(r)-t„)2 a^^dr 

J {x,tm) 
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and (3.6) we obtain 

(5.8) ||Ata^(t^+i) - [a^itm+i) -a^(tm;jiiL2 

"{^.'jtm+l) 




v/|x(r) + (t(r) -t„)2 a^^dr 



dx 



<At / Vlv^lWOP + l 



ai^dr 



< CAt"||v^(Wi)|Uo 




(a;,im) 
(a^,tm+i) 

(S,tm) 



where a^itm) '■= «^(a^,^m)- Since 

("^/i,) ^m) '-'^ ^m) 

then 

(5.9) \\aUtm)-a%trr.)\\l, 



/ i:)a^(x,, + s(x - a;/j),tm) ■ (x - x,,)(is, 
Jo 



At' 




Da%Xh + S{X - Xh),tm) ■ (V™ - \%trn))ds 



dx 



<Ce-^Af\\^r^-^r%t 



Using ( |5.8[ )-(5.9), we can bound the first term on the right-hand side of (5.7) as follows: 

(5.10) (At<(Wi) - [a%t^+i)~al{tm)],C''') 

< cAt^(ii«..|ii.([,„,„,,]xM3) + ^-1v- - v^(t„)lliO + Iwc^'Wh- 

To bound the second term on the right-hand side of ( 5.7[ ), writing 

a;"^+i(x)-^™(a;,) 

= (CJ"+1(X) - CU™(X)) + (CJ"(X) - UJ"\X)) + (CU'"(X) - Uj'''{Xh)), 

we then have 

(5.11) 11.;"^+^ - u-^Wh < At||^,|li^{[t„,wt]xM3). 
Next, it follows from 

a;™(x) - cu™(x) = At[ Du'^^x + s{x - x)) ■ ^'{t^)ds 

Jo 

that (set uJ"^ := uj'^{x)) 

(5.12) ||cu'"-cJ"||2, < CAt2||v^(t^)||i^||u;™||^i < CAt2||cu"^||^,. 
Finally, using the identity 

^"'{x) - uj'^ixh) = At[ Du^'ix + s{xh - x)) ■ (v=(t„) - v™)rfs 



20 



X. FENG AND M. NEILAN 



we get 
(5.13) 



2||, ,m||2 



|v^(U-v,"^lli. 



< CAt^\\^r'{t, 



,m\\2 



Combining (5.11 )-(5.13), we then bound the second term on the right-hand side of 
( 5.7[ ) as follows: 



(5.14) {u;-^'-uJ-,C^') < c(At||^,||i.([,^,^^^],^3) + Atl^™||^, 

+ At^K(U-v™||i.)+l||r^i||i.. 

Step 3: To get a lower bound of {^"'+^ - ^",^'"+^), let F^(x) := x - Atv;[^(a;). We 
then have 

det( J^„) = 1 + At^ (1 + i^ZA. - - (Cx. + C)) ' 

where Jp^ denotes the Jacobian of F^, and we have omitted the subscript h for notational 
convenience. Letting Ato = 0{6), we can conclude from the induction hypotheses that 
for At < Ato, Fm is invertible and det(J^-i) = 1 + Ce~'^At^. From this result we get 



(5.15) 
Thus, 

(5.16) 



linii. = (i + ^"^At^)liriii.. 



(iir+'iii^-(i + ^"'At')iiriiiO- 



and Remark 4.1 yield 



Step 4- Combining (5.7), (5.10), (5.14), (5.16), and using the induction hypotheses 



lie 



m+l||2 ||fm||2 
IIl2 ~ lis IIl2 



<Ce-^Mt^fll<Jli.([,^,„,,]xj 



:3) + 11^ \\m + 



\K-^'itm)\\h) 



+ At\\u;tf 



+At'\\n\h] 



< Ce 



-2 



{A.(| 



a 



e ||2 



\L'^{[tm.,tm + l]x1l 



-L 11/ ['"ll^ 

+ 11^ \\h^ 



+ CKe){C!{e)h''-'\\r{tJ\\l, + {Cl{e)h' + l)\\a%tj - a^U)) 



+ At|i^,iii.a*.,wi]xM3) + Atiriii.}. 



It follows from the inequality 



i«'(^m) - «niH- < Wit J - a^L'^ < iir iil2 + ii^'^iil^ 



that 



iir+ii^ 



lie 



m II 2 

Il2 



< 



Ce-'\ At'(\\aU\hat^,t^,,]><m + {Cl{e)K' + l)||a;' 



m||2 



+ C,^(.)C|(5)/^^^-lV^^(U||^. + At||c.,||i.([,^,^^^],, 



+ 



{ci{e)h' + i)Ae\\n\%]. 
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Applying the summation operator X]m=o noting that = we get 



2€-2|i ;,e||2 



\L^{[0,T];Hl) 



+ 



L^{[0,T];m) 



L2([0,T]xIR3) 



m\\2 

L2 (i 



m=0 



which by an apphcation of the discrete Gronwall inequahty yield 



(5.17) ||e'+iL2 < Ce^' (1 + E-\C,ie)h + l)At)'+' |At|| 



"rTlk2([0,T]xn 



+ VAt 



{C4{e)h + l)\\uj\\L2([o,T];m) 



\l2{[0,T];H'!) + ||^^t||L2([0,T]) 



We note hi = OjC^^ie )) = 0{el). Thus, for h < mm{hi, h2, h^} and At < 
min{Ato, /i^}, we have from (5.17), the triangle inequality, and (5.1) that 



(5.18) 



a%tk+i) - a/ < C5(£)|At||a^^||i2([o_T]xM3) + /i^[||a^||L2([o,T];//j) 

+ ll«t I|l2([o,T];HJ)] + C6{e)h^\\'ljj'^\\L2([o,T];m) (■ 



Thus, by Remark 4.1, we obtain the following estimates: 



(5.19) 



(5.20) 



\h2 < C7(£:)|At||a^^||L2([o_T]xIR3) + h^[\\a^\\L'2([0,T];H:>) 
+ I|l2([0,T];HJ)] + Ce{e)h^''^\\lp''\\L2(^[o^T];Hi)Y 

'(4+i) - i^f^^Wm < Cs{e)^^At\\aUL^[o,T]xR3) + y [\\a'\\L2{%T];W) 

+ Il«t I|l2([o,t];Hj)] + C&{e)h^~^\\'^p^\\L2^^^Q^'r]■,H^)^■ 
Step 5: We now verify the induction hypotheses. Set 

Cu{e) = C5{e)Ce{e)\\ij'\\L2(^[o,T];Hi), 

and let 



= Oi min 



CJe) 



Ati = O 



Ciiie 
mm{Cg{e),e'^h^} 

C5(£^)||"rTllL2([0,T]xM3) 



(Mel 



Cioie] 



2 

2j-3 



On noting that Ati < Ato, it follows from (5.18) that for h < min{/ii, /i2, h^, hr,} and 
At < min{Ati,/i2} 



a^(tfe+i) -a^+ii2 < Cmin{£^£^/l^C;^(£)}. 
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Thus, the first induction hypothesis (5.5) holds. 
Finally, let 



ha=0 



Ce{e)C7{e)\\ilj^\\L2(^lo^T];H^) 



2 

21-7 



by the definitions of /is, hf and Ati, (5.19), ( 3.6[) , and the inverse inequality we have for 
h < min{/ii, /i2, h^, h^, h^} and At < minjAti, h'^} 



L2 



< Ce 



-1 



h-^C7(e)lAt\\a 



|L2([0,T]xK3) 



< Ce 



Therefore, the second induction hypothesis (5.6) holds, and the proof is complete by 
setting /13 = min{/ii, /12, h^, h^}. □ 

Remark 5.1. In the two dimensional case, 



Ati = O 



C5(£^)||arrlU2([o,r]xM2) 



Remark 5.2. Recalling the definitions ofV^ and W^, we require k > r — 2 in order 
to obtain optimal order error estimate for tp^ in the H^-norm. 

6. Numerical experiments. In this section we shall present three 2-d numerical exper- 
iments. The first two experiments are done on the domain U = (0, 1)^, while the third 
experiment uses U = (0, 6)^. In all three experiments the fifth degree Argyris plate finite 
element (cf. [8]) is used to form V'\ and the cubic Lagrange element is employed to form 
W^. We recall that (see Section \2h the 2-d geostrophic fiow model has the exact same 
form as ( 1.1[ )-( |L5 ) except v and in ( |1.5 ) and (1.9) are replaced respectively by 



6.1. Test 1. The purpose of this test is twofold. First, we compute and ipj^ to view 
certain properties of these two functions. Specifically, we want to verify > and that 
'0/^ is strictly convex for m = 0, 1, M. Second, we calculate \\ip* — ipl\\ and ||a — 



for fixed h = 0.023 and At = 0.0005 in order to approximate \\^p* — ijj^ 



We set to solve problem (3.13)-(3.15) with the right-hand side of (3.15) being replaced 
by (F,Wh), and Vi and Wq being replaced by Vg^ and Wg^, respectively, where 



and II a — I 



^ = 9n, {vh, 1) = c(t) 

We use the following test functions and parameters 

gN{x,t) = te*(^?+^2)/2(a,^^^^,^ +a;2Z/:,.J, 
gD{x, t) = t\l + t{x\ + x2))e*("?+"2), 
F{x, t)=t{2 + At{xl + xl) + e{xl + x^)2)e*(^?+^i) 



c(t):=(V^*,l), 
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so that the exact solution of (1.1)-(1.5) is given by 



We record the computed solutions and plot the errors versus e in Figure 6.1 at tm = 
0.25. The figure shows that \\'ip*{tm) — V'^TIIj^^ — 0{ei), and since we have set both h 
and At very small, these results suggest that \\'>p*{tm) — 'ip'^{tm)\\m = 0{e^). Similarly, 

3 

we argue \\il)*{tm) - ^^(tm)||i/i = 0{ei) and \\i)*{tm) - V'^(im)||L2 = 0{e) based on our 
results. We note that these are the same convergence results found in [181 UHl EHl 125] . 
where the single Monge-Ampere equation was considered. We also notice that this test 
suggests that ||a(tM) — tth(^m)||L2 may not converge, which suggests that the convergence 
can only be possible in a weaker norm such as . 

Next, we plot a"^, and ip'f^ for = 0.1, 0.4, and 1.0 with h = 0.05, At = 0.1 in Figure 



6.2 The figure shows that a^(tm) > and the computed solution ifj'^ is clearly convex 
for all tr. 




Fig. 6.1. Test 1: Change of \\ip' (tm) - '/'"ll w.r.t. e. h = 0.023, At = 0.0005, = 0.25. 



and 



6.2. Test 2. The goal of this test is to calculate the rate of convergence of H^"^ — 4'', 



for a fixed e while varying At and h with the relation At = h 



(3.13)-(3.15) but with a new boundary condition: 



du 



Let V;'^ and W^^ 



. We solve 
be defined 
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Test 1: Computed t_m=0.1 Max: 1.105 Test 1: Computed a. t_m=0.1 Max: 0.0170 




Min: 0.999 Min: 0.0100 

Test 1: Computed qj'', t_m=0.4 Max: 1.499 Test 1: Computed a, t_m:^0.4 Max: 0.641 




Min: 0.995 Min: 0.160 

Test 1: Computed qj*, t_m=1.0 Max: 2.731 Test 1: Computed a, t_m = 1.0 Max: 22.157 




Fig. 6.2. Computed -0^ (^^fi) and (right) for Test 1 at tm — 0.1 (top), tm — 0.4 (middle), and tm — 1-0 (bottom). 
At ^ 0.1, h ^ 0.05. 

in the same way as in Test 1 using the following test functions and parameters 

c(^) = (v^^l), 

= t2(i + t(a;2 + x2))e*("?+"2)^ 
F = t{2 + At{xl + xl) + t\xl + a;2)2)e*(-'?+"i) 

- ^p*(-?+-l)/2/Q2 + 56r.T? + xht + 16^2. 2 , ^2.2 , ,3/2 , 2^3^ 
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h 


At 


\\r{trr.)-^r\\L^ 








0.08333 


0.00694 


0.000214135 


0.000978608 


0.004434963 


0.003456864 


0.05 


0.0025 


6.15715E-05 


0.000281367 


0.001274611 


0.001009269 


0.03066 


0.00094 


1.42185E-05 


6.49825E-05 


0.000294575 


0.000232896 


0.02384 


0.00057 


7.13357E-06 


3.25959E-05 


0.000147586 


0.000116731 



Table 1: Change of ||^A=(tM) - V," || w.r.t. At = . e = 0.01, = 0.25. 



SO that the exact solution of ( 1.6 )-( 1.12) is given by 

a%x, t) = t\l + t{xl + x2))e*("?+"') - rf2gt(.-?+-i)/2(8 + 8t(x? + xl) + t^xj + xj)^). 

The errors at time = 0.25 are hsted in Table 1 and are plotted verses At in Figure 
6.3 The results clearly indicate that ||a'^(trrt) — ct/TlU^ = 0{At) and \\ip'^itm) — i'h'W = 



0{At) in all norms as expected by the analysis in the previous section. 







Fig. 6.3. Test 2: Change of il-0^(tAf) - ipff \\ w.r.t. At = /i^ . e = 0.01, t„, = 0.25. 



6.3. Test 3. For this test, we solve problem (3.13)-(3.15) with domain U = (0, 6)^ and 
initial condition 



«o(a;) = 3 X[2,4]x [2.25,3.75] (4 - Xi){xi - 2)(3.75 - X2)ix2 - 2.25), 
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where X[2,4]x [2.25,3.75] denotes the characteristic function of the set [2,4] x [2.25,3.75]. We 
comment that the exact solution of this problem is unknown. We plot the computed 
and ip'^ at times = 0, = 0.05, and tm = 0.1, and tm = 0.15 in Figure 6.4| with 
parameters At = 0.001, h = 0.05, and e = 0.01. As expected, the figure shows that 
> and is convex for all m. 



Test 3: Computed a, t_m==0.0 



Test 3: Computed t m^O.O 




Test 3: Computed a, t_m-0. 05 



Test 3: Computed «, t_m=0.1 



Min: -1.77e-3 
Max: 0.0695 





Min: -7.382e-3 
Max: 0.0690 

0.06 




Test 3: Computed ly, t_m=0.1 




Min: -7.704e-3 



Max: 0.0522 

10.05 
... 



-0.01 


1-0.01 

-0.02 

-0.03 
Min: -0.0332 



Fig. 6.4. Test 3: Computed ol^ (left) and (right) at tm — (top), tm — 0.05 (middle), and tm — 0.1 (bottom). 

At = 0.01, h = 0.05, £ = 0.01 
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